{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "74F0A5FB029B49568E19048977C348C3",
        "id": "0187CA1FFA87432784F6F2DF788F9F42",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "# <center> Lecture11 : Evaluating Regression Models </center>  \n",
        " \n",
        "## <center> Instructor: Dr. Hu Chuan-Peng </center> "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "5BF0DA15FE764C37B8B6785557411E64",
        "id": "539711953E154C04A3B560AC42FA1E36",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "在Lecture 10的末尾，我们引入了是否吸烟这个变量，并建立了一个新的回归模型。  \n",
        "\n",
        "在本节课中，我们将讨论更复杂的模型，例如同时包含两个预测变量(压力和是否吸烟)，甚至是两个变量的交互项。  \n",
        "\n",
        "🤔这节课需要关注的问题：  \n",
        "1. 如何评估不同的模型？哪个模型最好？  \n",
        "2. 什么时候该使用什么模型？  \n",
        "3. 贝叶斯分析流程中各部分的意义。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "D238022A3BC24276B0C7FE2BBDF95D34",
        "id": "55C10C12AB6E4F92A88A8E594B3BBE8F",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 模型回顾  \n",
        "\n",
        "在之前的课程中，我们分别以压力和是否吸烟作为预测因子定了一般线性模型  \n",
        "\n",
        "$$  \n",
        "Y_i | \\beta_0, \\beta_1, \\sigma \\stackrel{ind}{\\sim} N\\left(\\mu_i, \\sigma^2\\right) \\;\\; \\text{ with } \\;\\; \\mu_i = \\beta_0 + \\beta_1X_i  \n",
        "$$  \n",
        "\n",
        "- 两种模型可以用相同的回归公式表示。  \n",
        "- 区别在于：  \n",
        "\t- 两者的先验不同。  \n",
        "\t- 对于离散变量中回归系数的解释不同。  \n",
        "\n",
        "<table>  \n",
        "        <tr>  \n",
        "\t\t\t\t\t<td>模型</td>  \n",
        "           <td>model1</td>  \n",
        "\t\t\t\t\t <td> model2</td>  \n",
        "        </tr>  \n",
        "        <tr>  \n",
        "            <td>自变量</td>  \n",
        "\t\t\t\t\t<td>压力(连续变量)</td>  \n",
        "\t\t\t\t\t<td>吸烟(离散变量)</td>  \n",
        "        </tr>  \n",
        "\t  <tr>  \n",
        "            <td>自变量含义</td>  \n",
        "\t\t\t\t\t<td>`0` 表示不吸烟，`1`表示吸烟</td>  \n",
        "\t\t\t\t\t<td>14-70的压力评级</td>  \n",
        "   </tr>  \n",
        "\t <tr>  \n",
        "            <td>先验</td>  \n",
        "\t\t\t\t\t<td>  \n",
        "β0 ~ N(45, 10) <br>  \n",
        "β1 ~ N(0, 10)  <br>  \n",
        "σ ~ Exp(0.6) <br>  \n",
        "\t\t </td>  \n",
        "\t\t\t\t\t<td>  \n",
        "β0 ~ N(45, 10)  <br>  \n",
        "β1 ~ N(0, 10)  <br>  \n",
        "σ ~ Exp(0.6) <br>  \n",
        "\t\t </td>  \n",
        "  </tr>  \n",
        "\t  <tr>  \n",
        "            <td>可视化</td>  \n",
        "\t\t\t\t\t<td><img src=\"https://cdn.kesci.com/upload/s3y4rmsym1.png?imageView2/0/w/600/h/600\" alt=\"\" width=\"100\" height=\"200\"> </td>  \n",
        "\t\t\t<td><img src=\"https://cdn.kesci.com/upload/s40nk3kagu.png?imageView2/0/w/600/h/600\" alt=\"\" width=\"100\" height=\"200\"> </td>  \n",
        "   </tr>  \n",
        "</table>  \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "_id": "98C3B096B8B842A48E261C2361BAD7AE",
        "collapsed": false,
        "id": "D9C003C53DD940828E5BE3C8D6892627",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 导入 pymc 模型包，和 arviz 等分析工具 \n",
        "import pymc as pm\n",
        "import arviz as az\n",
        "import seaborn as sns\n",
        "import scipy.stats as st\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import xarray as xr\n",
        "import pandas as pd\n",
        "\n",
        "# 忽略不必要的警告\n",
        "import warnings\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "_id": "E7221776DE4344B0983650BBC8DF3E5A",
        "collapsed": false,
        "id": "6BC6C8DDD31444858622BD03349FD953",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 通过 pd.read_csv 加载数据 Data_Sum_HPP_Multi_Site_Share.csv\n",
        "try:\n",
        "  df_raw = pd.read_csv('/home/mw/input/bayes20238001/Data_Sum_HPP_Multi_Site_Share.csv')\n",
        "except:\n",
        "  df_raw = pd.read_csv('data/Data_Sum_HPP_Multi_Site_Share.csv')\n",
        "\n",
        "# 筛选站点为\"Tsinghua\"的数据\n",
        "df = df_raw[df_raw[\"Site\"] == \"Tsinghua\"]\n",
        "\n",
        "# 选择需要的两列\n",
        "df = df[[\"stress\",\"scontrol\",\"smoke\"]]\n",
        "\n",
        "#1 表示吸烟，2表示不吸烟\n",
        "df[\"smoke\"] =  np.where(df['smoke'] == 2, 0, 1)\n",
        "df[\"smoke_recode\"] =  np.where(df['smoke'] == 1, \"yes\", \"no\")\n",
        "\n",
        "#设置索引\n",
        "df[\"index\"] = range(len(df))\n",
        "df = df.set_index(\"index\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "AD2625CBC10D460CA4DF8C8845234E1E",
        "id": "6E908DCA7A6E44EBB4D48ECD484D49FB",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**快速运行模型定义、后验参数采样、后验预测分布部分**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "_id": "1CBE9B76A99C46FF808B0D02CE75FF2D",
        "collapsed": false,
        "id": "EF42B25098AC4080978B6AF5326A5C5B",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model() as model1:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=45, sigma=10)        #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)         #定义beta_1\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                 #定义sigma\n",
        "\n",
        "    x = pm.MutableData(\"smoke\",df.stress, dims=\"obs_id\")                     #x是自变量压力水平\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + beta_1*x, dims=\"obs_id\")            #定义mu，讲自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol, dims=\"obs_id\")\n",
        "\n",
        "\n",
        "with pm.Model() as model2:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=45, sigma=10)        #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)         #定义beta_1\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                 #定义sigma\n",
        "\n",
        "    smoke = pm.MutableData(\"smoke\",df.smoke, dims=\"obs_id\")            #x为是否吸烟，0 1编码。\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + beta_1*smoke, dims=\"obs_id\")  #定义mu，讲自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol, dims=\"obs_id\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "_id": "893EDF1EAE0944A7811C8CA848079F22",
        "collapsed": false,
        "id": "AC411D63FA6F444892725791E3AD5101",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "pm.model_to_graphviz(model1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "_id": "01E1F363D57E435D8A8AE5154ECE5DE1",
        "collapsed": false,
        "id": "9C58AE6270B542D8AF813FDE8972B77A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "#===========================\n",
        "#     注意！！！以下代码可能需要运行5分钟左右\n",
        "#===========================\n",
        "with model1:\n",
        "    model1_trace = pm.sample(draws=5000,                  # 使用mcmc方法进行采样，draws为采样次数\n",
        "                      tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                      chains=4,                     # 链数\n",
        "                      discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                      random_seed=84735)\n",
        "\n",
        "with model2:\n",
        "    model2_trace = pm.sample(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "                        tune=1000,                   # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                        chains=4,                    # 链数\n",
        "                        discard_tuned_samples= True, # tune的结果将在采样结束后被丢弃\n",
        "                        random_seed=84735)                     "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "_id": "EB054F6D1CDE4F88910D1C0965439041",
        "collapsed": false,
        "id": "B5041DF6EA33432E91792241E1BBC431",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 后验预测\n",
        "with model1:\n",
        "    model1_ppc = pm.sample_posterior_predictive(model1_trace, random_seed=84735)\n",
        "with model2:\n",
        "    model2_ppc = pm.sample_posterior_predictive(model2_trace, random_seed=84735)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "00A86C01F50146179D3895B0D8D61873",
        "id": "549EE74A94F94C0EB96C256ACBCC0A33",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 评估后验预测  \n",
        "\n",
        "这里我们仍然使用上节课介绍过的两个指标：**MAE**与**后验预测区间**。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "1E1E8DBA4E2046DA8215AFF2AC7A4CF8",
        "id": "C2FC7565FB8A43068E37685B46264876",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### MAE  \n",
        "\n",
        "$$  \n",
        "\\text{MAE} = median\\left(|Y_i - \\hat{Y}_i|\\right)  \n",
        "$$  \n",
        "\n",
        "MAE的结果显示：  \n",
        "- 模型1的预测误差小于模型2，即模型1的预测能力优于模型2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "_id": "B8688466D8004E6A85AA5555348180FB",
        "collapsed": false,
        "id": "089ED3C29B1747269C8FB38AA5FA13D2",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义计算 MAE 函数\n",
        "from statistics import median\n",
        "def MAE(model_ppc):\n",
        "    # 计算每个X取值下对应的后验预测模型的均值\n",
        "    pre_x = model_ppc.posterior_predictive[\"y_est\"].stack(sample=(\"chain\", \"draw\"))\n",
        "    pre_y_mean = pre_x.mean(axis=1).values\n",
        "\n",
        "    # 提取观测值Y，提取对应Y值下的后验预测模型的均值\n",
        "    MAE = pd.DataFrame({\n",
        "        \"scontrol_ppc_mean\": pre_y_mean,\n",
        "        \"scontrol_original\": df.scontrol\n",
        "    })\n",
        "\n",
        "    # 计算预测误差\n",
        "    MAE[\"pre_error\"] = abs(MAE[\"scontrol_original\"] -\\\n",
        "                            MAE[\"scontrol_ppc_mean\"])\n",
        "\n",
        "    # 最后，计算预测误差的中位数\n",
        "    MAE = median(MAE.pre_error)\n",
        "    return MAE\n",
        "\n",
        "# 输出结果\n",
        "model1_MAE = MAE(model1_ppc)\n",
        "print(f\"模型1(连续变量)MAE: {model1_MAE:.2f}\")\n",
        "model2_MAE = MAE(model2_ppc)\n",
        "print(f\"模型2(离散变量)MAE: {model2_MAE:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "CBED90DC4FE14094A8D81872F16F1818",
        "id": "FCC68730726741669D39838681E9F0F4",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 后验预测区间  \n",
        "\n",
        "结果发现，两个模型中超出95%hdi的观测值的数量为9。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "_id": "899D635CD2264846886CA9D30522A056",
        "collapsed": false,
        "id": "51BB9F321F3240BE8EBCBF10D0D4E72B",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义\n",
        "def counter_outlier(model_ppc, hdi_prob=0.95):\n",
        "    # 将az.summary生成的结果存到hdi_multi这个变量中，该变量为数据框\n",
        "    hdi = az.summary(model_ppc, kind=\"stats\", hdi_prob=hdi_prob)\n",
        "    lower = hdi.iloc[:,2].values\n",
        "    upper = hdi.iloc[:,3].values\n",
        "\n",
        "    # 将原数据中的自我控制分数合并，便于后续进行判断\n",
        "    y_obs = model_ppc.observed_data[\"y_est\"].values\n",
        "\n",
        "    # 判断原数据中的压力分数是否在后验预测的95%可信区间内，并计数\n",
        "    hdi[\"verify\"] = (y_obs <= lower) | (y_obs >= upper)\n",
        "    hdi[\"y_obs\"] = y_obs\n",
        "    hdi_num = sum(hdi[\"verify\"])\n",
        "\n",
        "    return hdi_num\n",
        "\n",
        "# 输出结果\n",
        "model1_outliers = counter_outlier(model1_ppc)\n",
        "print(f\"模型1(连续变量)超出95%hdi: {model1_outliers:.2f}\")\n",
        "model2_outliers = counter_outlier(model2_ppc)\n",
        "print(f\"模型2(离散变量)超出95%hdi: {model2_outliers:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "09D64CC54CDD496E93425CF35C18D782",
        "id": "3E3AEE345B70459796FD76CC79064274",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 后验预测分布  \n",
        "\n",
        "* 通过 arviz 绘制两个模型的后验预测分布  \n",
        "* 可见，两个模型对于数据总体的预测性能相当"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {
        "_id": "03D90670E2874BF5A24187B2A332860E",
        "collapsed": false,
        "id": "60638F7CC254478CA398E8312531E583",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义两个模型ppc可视化的画布\n",
        "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(20,10))\n",
        "\n",
        "# 绘制模型1 ppc\n",
        "az.plot_ppc(model1_ppc, num_pp_samples=500, ax=axes[0])\n",
        "axes[0].set_title(\"model1 (stress)\")\n",
        "# 绘制模型2 ppc\n",
        "az.plot_ppc(model2_ppc, num_pp_samples=500, ax=axes[1])\n",
        "axes[1].set_title(\"model2 (smoke)\")\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "59F55A34294640119090B2CD91EDA820",
        "id": "878F5F69394341388D340773D2963CFD",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 多元线性回归：增加预测变量  \n",
        "\n",
        "之前的两个模型，我们分辨考虑两个预测变量对于模型的预测性能。  \n",
        "\n",
        "🤔 是否可以将吸烟情况和压力同时加入到模型中?  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "A932777D911B4CA2970DD761F950B222",
        "id": "195E468E885B4A18B404B8FECF51A61A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型定义  \n",
        "\n",
        "$$  \n",
        " \n",
        "\\begin{array}{lcrl}  \n",
        "\n",
        "\\text{data:} & \\hspace{.05in} & Y_i|\\beta_0,\\beta_1,\\beta_2,\\sigma & \\stackrel{ind}{\\sim} N(\\mu_i, \\; \\sigma^2)  \\;\\; \\text{ with } \\;\\; \\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2} \\\\  \n",
        "\\text{priors:} & & \\beta_{0}  & \\sim N\\left(50, 10^2 \\right)  \\\\  \n",
        "                    & & \\beta_1  & \\sim N\\left(0, 10^2 \\right) \\\\  \n",
        "                    & & \\beta_2  & \\sim N\\left(0, 10^2 \\right) \\\\  \n",
        "                    & & \\sigma   & \\sim \\text{Exp}(0.6)  .\\\\  \n",
        "\\end{array}  \n",
        "\n",
        "$$  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "CDAC3964400F4FB5A40FD87949F6D85C",
        "id": "5F275FEC0FF346F29EC1FB798863D84D",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**各参数(回归系数)的意义：**  \n",
        "\n",
        "* 在自变量中，$X_{i2}$为离散变量，`0` 表示不吸烟，`1`表示吸烟  \n",
        "\n",
        "$$  \n",
        "X_{i2} = \\begin{cases}  \n",
        "1 & \\text{Yes} \\\\  \n",
        "0 & \\text{No} \\\\  \n",
        "\\end{cases}  \n",
        "$$  \n",
        "\n",
        "在之前自变量仅为离散变量的例子中：  \n",
        "\n",
        "* 当$X_{i2} = 0$时，$\\beta_0  =  \\beta_0 + \\beta_1 \\cdot 0$  \n",
        "\n",
        "    *  $\\beta_0$ 表示不抽烟情况下，自我控制分数的均值  \n",
        "\n",
        "* 当$X_{i2} = 1$时 $\\beta_0 + \\beta_1 \\cdot 1 = \\beta_0 + \\beta_1$  \n",
        "\n",
        "    * $\\beta_0 + \\beta_1$ 表示抽烟情况下，自我控制分数的均值；则$\\beta_1$表示了从不抽烟到抽烟，自我控制分数的变化程度  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "14EA7AA2C32742819311263026A26CCF",
        "id": "5DA310298F6A468C87802DC5CF737935",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**各参数(回归系数)的意义：**  \n",
        "\n",
        "在现在的例子中，自变量包括离散变量和连续变量。  \n",
        "\n",
        "* 当$X_{i2} = 0$时，$\\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2 \\cdot 0 = \\beta_0 + \\beta_1X_{i1}$  \n",
        "\n",
        "    * 表示不吸烟情况下，自我控制分数随压力分数变化的情况，二者的关系可以被简化为一条直线。  \n",
        "\n",
        "* 当$X_{i2} = 1$时，$\\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2 \\cdot 1 = (\\beta_0 + \\beta_2) + \\beta_1X_{i1}$  \n",
        "\n",
        "    * 表示吸烟情况下，自我控制分数随压力分数变化的情况。  \n",
        "\n",
        "> 注意，在该模型中，我们假设在两种吸烟条件下，压力对自我控制的影响是相同的，即不存在吸烟对于压力的调节作用。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "C79007F316584A2C83555CF4975BC130",
        "id": "61B741F0FB37404AA5B2FFC0E3DEF361",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "🤔思考：  \n",
        "\n",
        "$\\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2}$  \n",
        "\n",
        "- 如果离散变量 $X_{i2}$改变编码，`-1` 表示不吸烟，`1`表示吸烟。此时各回归系数的含义是什么？  \n",
        "\n",
        "![](https://pic3.zhimg.com/80/v2-f95e51b7e6e26feab9ca91077d9c244e_1440w.webp)  \n",
        "\n",
        "> 延展阅读：  \n",
        "> 线性模型中无序因子变量的对比方式与回归系数的关系：原理，困境及其解决 https://zhuanlan.zhihu.com/p/88270555  \n",
        "> 通过线性模型详解方差分析中“对比”（Contrasts）的数理原理 https://zhuanlan.zhihu.com/p/119394473"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "13232A80CDC34CAAB3F35CE9DCBE2EAF",
        "id": "48E6BA887B4042D59431A9D670F0E96A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**直接运行模型定义、MCMC采样、后验预测部分**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "_id": "CF0D7F4E507742AA9F5884B7D223867E",
        "collapsed": false,
        "id": "0D94AF6C5EAC4D4E9857D5A5B6F261B1",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model() as model3:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=50, sigma=10)          #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)           #定义beta_1\n",
        "    beta_2 = pm.Normal(\"beta_2\", mu=0, sigma=10)           #定义beta_2\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                   #定义sigma\n",
        "\n",
        "    stress = pm.MutableData(\"stress\",df.stress, dims=\"obs_id\")                           #自变量1:压力水平\n",
        "    smoke = pm.MutableData(\"smoke\",df.smoke, dims=\"obs_id\")                              #自变量2:是否吸烟\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + beta_1*stress + beta_2*smoke, dims=\"obs_id\")    #定义mu，自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol, dims=\"obs_id\") "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "_id": "C8A527EEBC534EDE8B0BF33BBED6B429",
        "collapsed": false,
        "id": "4FB12A12C650442986A46754A8B4E1CA",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "#===========================\n",
        "#     注意！！！以下代码可能需要运行1-2分钟左右\n",
        "#===========================\n",
        "with model3:\n",
        "    model3_trace = pm.sample(draws=5000,            # 使用mcmc方法进行采样，draws为采样次数\n",
        "                      tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                      chains=4,                     # 链数\n",
        "                      discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                      random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "_id": "E8926B763DD04A4DAA79386771C3C17A",
        "collapsed": false,
        "id": "B41C37A9609B4E71BCECF9FEA6825FB7",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with model3:\n",
        "    model3_ppc = pm.sample_posterior_predictive(model3_trace, random_seed=84735) "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "E887D58D309543CA9C0ACFF5E32F6092",
        "id": "CC55CA01672B4193B539F60F010153CA",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 可视化后验预测结果  \n",
        "\n",
        "#### 后验预测的线性模型  \n",
        "\n",
        "我们可以分别查看**吸烟情况**与**压力**对**自我控制**的预测情况  \n",
        "\n",
        "1. $\\mu_i = \\beta_0 + \\beta_1X_{i1}$：自我控制分数随压力分数变化的情况  \n",
        "\n",
        "2. $\\mu_i = \\beta_0 + \\beta_2X_{i2}$：自我控制分数随吸烟变化的情况  \n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "_id": "01A28E9EEB5B412998BBD100F3324654",
        "collapsed": false,
        "id": "7D85C5423653492B908C6E8975E6883A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.summary(model3_trace, kind = \"stats\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "_id": "CE5019B381664D88874C7D8C70D31C0B",
        "collapsed": false,
        "id": "E35BAAB17ABB4CCC87AAA857F50F1B38",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 将参数beta_0, beta_1 与自变量 压力分数 进行组合，生成回归线\n",
        "model3_trace.posterior[\"stress_ymodel\"] = model3_trace.posterior[\"beta_0\"] +\\\n",
        "                                              model3_trace.posterior[\"beta_1\"] * model3_trace.constant_data[\"stress\"]\n",
        "\n",
        "# 将参数beta_0, beta_2 与自变量 吸烟情况 进行组合，生成回归线\n",
        "model3_trace.posterior[\"smoke_ymodel\"] = model3_trace.posterior[\"beta_0\"] +\\\n",
        "                                              model3_trace.posterior[\"beta_2\"] * model3_trace.constant_data[\"smoke\"]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "metadata": {
        "_id": "41C2C1A3D3094A83BD131C40DA4C61C2",
        "collapsed": false,
        "id": "89FFE43C0E7940308BDEE375CEA3310D",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义两幅图的画布\n",
        "fig, axes = plt.subplots(1,2,figsize=(18, 8))\n",
        "\n",
        "#----------------------------------------#\n",
        "#绘制预测变量为压力情况下的后验预测线性模型\n",
        "az.plot_lm(idata=model3_trace,\n",
        "           y= \"y_est\",\n",
        "           x= \"stress\",\n",
        "           y_model = model3_trace.posterior[\"stress_ymodel\"],\n",
        "           kind_model = \"hdi\",\n",
        "           textsize=15,\n",
        "           grid=False,\n",
        "           axes=axes[0],\n",
        "           legend=False,\n",
        "           y_model_fill_kwargs = {\"color\":\"grey\"})\n",
        "#设置坐标轴标题\n",
        "axes[0].set_xlabel('Stress', fontsize=16) \n",
        "axes[0].set_ylabel('Self control', fontsize=16) \n",
        "#----------------------------------------#\n",
        "#绘制预测变量为吸烟情况下的后验预测线性模型\n",
        "az.plot_lm(idata=model3_trace,\n",
        "           y= \"y_est\",\n",
        "           x= \"smoke\",\n",
        "           y_model = model3_trace.posterior[\"smoke_ymodel\"],\n",
        "           kind_model = \"hdi\",\n",
        "           textsize=15,\n",
        "           grid=False,\n",
        "           axes=axes[1],\n",
        "           y_model_fill_kwargs = {\"color\":\"grey\"})\n",
        "\n",
        "#设置坐标轴标题\n",
        "axes[1].set_xlabel('Smoke', fontsize=16) \n",
        "axes[1].set_ylabel('Self control', fontsize=16) \n",
        "#----------------------------------------#\n",
        "# 设置坐标轴字体大小、图例位置、总标题\n",
        "axes[0].tick_params(labelsize=14)  \n",
        "axes[1].tick_params(labelsize=14)\n",
        "axes[1].legend(loc=\"upper right\", fontsize=12)\n",
        "\n",
        "fig.suptitle('Posterior predictive regression models', fontsize=20)\n",
        "sns.despine()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "D6427797400648F59665506F02BF5DED",
        "id": "6BBA6656248C479AAE43114A45F99039",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "结果发现：  \n",
        "* 对于吸烟情况的预测出现了偏离，原因在于加入变量压力会影响对于吸烟效应的估计。  \n",
        "* 可以发现，加入压力变量后，吸烟的效应受到了一定的影响 (-1.75 vs. -1.06)。  \n",
        "* 但更主要的是，模型截距的变化很大 (55.82 vs. 39.55)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "_id": "2A1AB823F18E479CA7C62F42897FBCCF",
        "collapsed": false,
        "id": "23208D521B4C4A9D94570AB059A3277D",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "result1 = az.summary(model2_trace, kind=\"stats\")\n",
        "result2 = az.summary(model3_trace, kind=\"stats\")\n",
        "print(f\"模型2中，截距:{result1.loc['beta_0','mean']:.2f}, 吸烟的回归系数:{result1.loc['beta_1','mean']:.2f}\")\n",
        "print(f\"模型3中，截距:{result2.loc['beta_0','mean']:.2f}, 吸烟的回归系数:{result2.loc['beta_2','mean']:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "30595F66B9E049B984B4D34CB20A31C3",
        "id": "6456EC139B2A4CAA8D320CDB9A998A9A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 绘制在不同吸烟状况下，压力对自我控制影响的回归模型  \n",
        "\n",
        "注意：该模型假设在两种吸烟条件下压力的预测能力相同，因此，两条回归预测线(深绿色)是平行的。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "_id": "3255FCD4199A4BA8989A89FFA55D6A95",
        "collapsed": false,
        "id": "ABA52F80567A4E6EBAA0DDDEADE84BBF",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 分别获取smoke=0及smoke=1时，压力变量对应的索引值\n",
        "no_smoke_index = xr.DataArray(df.query(\"smoke == 0\")[\"stress\"])[\"index\"].values\n",
        "smoke_index = xr.DataArray(df.query(\"smoke == 1\")[\"stress\"])[\"index\"].values"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "_id": "8F94560D83D14659A1740CC3B71CB55D",
        "collapsed": false,
        "id": "9636EBCB24AE4DC18E70036B9563A4DB",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 根据上一步获得的索引值，将参数的后验预测分布分成两部分\n",
        "# 一部分为吸烟情况下，参数的后验预测分布、压力分数、自我控制分数\n",
        "# 另一部分为不吸烟情况下，参数的后验预测分布、压力分数、自我控制分数\n",
        "nosmoke_para = model3_trace.sel(obs_id=no_smoke_index)\n",
        "smoke_para = model3_trace.sel(obs_id=smoke_index)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "_id": "9085D9854A1C4C30A233D516C377C39C",
        "collapsed": false,
        "id": "E065C508677A411B9DA6F32031DAFEBE",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 在不吸烟情况下，将参数与自变量进行组合，生成回归线\n",
        "nosmoke_para.posterior[\"y_model\"] = nosmoke_para.posterior[\"beta_0\"] +\\\n",
        "                                    nosmoke_para.posterior[\"beta_1\"] * nosmoke_para.constant_data[\"stress\"] +\\\n",
        "                                    nosmoke_para.posterior[\"beta_2\"] * nosmoke_para.constant_data[\"smoke\"]\n",
        "\n",
        "# 在吸烟情况下，将参数与自变量进行组合，生成回归线\n",
        "smoke_para.posterior[\"y_model\"] = smoke_para.posterior[\"beta_0\"] +\\\n",
        "                                  smoke_para.posterior[\"beta_1\"] * smoke_para.constant_data[\"stress\"] +\\\n",
        "                                  smoke_para.posterior[\"beta_2\"] * smoke_para.constant_data[\"smoke\"]\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "_id": "5CE0FF25D02449DA930A6D184346AAA5",
        "collapsed": false,
        "id": "DE85D43A933642AE9C809D55C4C94024",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "def plot_ppc_lm(nosmoke_para, smoke_para):\n",
        "    fig, axes = plt.subplots(figsize=(10, 8))\n",
        "\n",
        "    #绘制smoke=0情况下的后验预测线性模型\n",
        "    az.plot_lm(idata=nosmoke_para,\n",
        "            y= \"y_est\",\n",
        "            x= \"stress\",\n",
        "            y_model = nosmoke_para.posterior[\"y_model\"],\n",
        "            kind_model = \"hdi\",\n",
        "            textsize=15,\n",
        "            grid=False,\n",
        "            axes=axes,\n",
        "            legend=False,\n",
        "            y_model_fill_kwargs = {\"color\":\"orange\"},         # 设置 hdi 区间颜色\n",
        "            y_kwargs =  {\"color\":\"orangered\",\"alpha\":0.3}     # 设置观测值颜色\n",
        "            )\n",
        "\n",
        "    #----------------------------------------#\n",
        "    #绘制smoke=1情况下的后验预测线性模型\n",
        "    az.plot_lm(idata=smoke_para,\n",
        "            y= \"y_est\",\n",
        "            x= \"stress\",\n",
        "            y_model = smoke_para.posterior[\"y_model\"],\n",
        "            kind_model = \"hdi\",\n",
        "            textsize=15,\n",
        "            grid=False,\n",
        "            axes=axes,\n",
        "            y_model_fill_kwargs = {\"color\":\"blue\"},         # 设置 hdi 区间颜色\n",
        "            y_kwargs =  {\"color\":\"navy\"},                   # 设置观测值颜色\n",
        "            legend=False)\n",
        "\n",
        "    #设置坐标轴标题\n",
        "    axes.set_xlabel('Stress', fontsize=16) \n",
        "    axes.set_ylabel('Self control', fontsize=16) \n",
        "\n",
        "    #设置图例\n",
        "    handles = [\n",
        "            plt.Line2D([], [], color=\"orange\",\n",
        "                        label = 'smoke = no'),\n",
        "            plt.Line2D([], [], color=\"blue\",\n",
        "                        label = 'smoke = yes')]\n",
        "\n",
        "    axes.legend(handles=handles)\n",
        "    #----------------------------------------#\n",
        "    # 设置坐标轴字体大小、总标题\n",
        "    axes.tick_params(labelsize=14)  \n",
        "    fig.suptitle('Posterior predictive regression models', fontsize=20)\n",
        "    sns.despine()\n",
        "\n",
        "plot_ppc_lm(nosmoke_para, smoke_para)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "561275E408AD4388992D2658C7AD1C53",
        "id": "22DFA6488B294BE19A5854CA4ADED096",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 后验预测分布"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "_id": "6404110E644745459A1A05BC14BEFC7B",
        "collapsed": false,
        "id": "950C6560EE5D431AA67FADCE78A045B2",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.plot_ppc(model3_ppc, num_pp_samples = 500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "1EA79453A64A4D1B9E9B2349DE8367D9",
        "id": "ADFACE35C9BA49799D51A2AD53E5AA3B",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 评估后验预测  \n",
        "\n",
        "#### MAE  \n",
        "\n",
        "$$  \n",
        "\\text{MAE} = median\\left(|Y_i - \\hat{Y}_i|\\right)  \n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "_id": "1BA1409F08ED43A4B38E16C33BFAF02B",
        "collapsed": false,
        "id": "06CBEA71E47B421E93EFB1053762C722",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "model3_MAE = MAE(model3_ppc)\n",
        "print(f\"模型3(多元回归)MAE: {model3_MAE:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "931D74594DAE4AC69BCE15E4503AFE6C",
        "id": "49FB07795413488CAA474D582D21EA64",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 后验预测区间  \n",
        "\n",
        "* 使用`az.summary`，可以直接得到每个Y后验预测值95%的可信区间  \n",
        "* 将真实数据与其合并，判断其是否在对应的95%可信区间内，并进行绘图  \n",
        "\n",
        "> 注意：在真实数据中，一个X可能对应着多个Y，因此可信区间在绘制时会出现重叠"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "_id": "6F2865B0ED8B4ED199152D2FF2AA54C2",
        "collapsed": false,
        "id": "F1CA38E6288B41E2A76FB37AD420011B",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "model3_outliers = counter_outlier(model3_ppc)\n",
        "print(f\"模型3(多元回归)超出95%hdi: {model3_outliers:.2f}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "_id": "45ECD112A399404A9E59740B5C838A1A",
        "collapsed": false,
        "id": "196234E14D484E4B8AF54D5820699F3E",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 定义绘图函数\n",
        "def plot_ppi(model_ppc):\n",
        "    fig, ax =  plt.subplots(figsize=(8,6))\n",
        "\n",
        "    # 将az.summary生成的结果存到hdi_multi这个变量中，该变量为数据框\n",
        "    hdi_multi = az.summary(model_ppc, hdi_prob=0.95)\n",
        "    hdi_multi = hdi_multi.reset_index()\n",
        "    # 将原数据中的压力分数与自我控制分数合并，便于后续进行判断\n",
        "    hdi_multi[\"x_obs\"] = model_ppc.constant_data[\"stress\"].values\n",
        "    hdi_multi[\"y_obs\"] = model_ppc.observed_data[\"y_est\"].values\n",
        "\n",
        "    # 绘制95%的可信区间\n",
        "    HDI = ax.vlines(hdi_multi[\"x_obs\"], \n",
        "            hdi_multi[\"hdi_2.5%\"], hdi_multi[\"hdi_97.5%\"], \n",
        "            color=\"orange\", \n",
        "            alpha=0.5,\n",
        "            label=\"94% HDI\")\n",
        "\n",
        "    #绘制真实值的散点图，并根据是否落在区间内设置观测值的不同颜色\n",
        "    colors = np.where((hdi_multi[\"y_obs\"] >= hdi_multi[\"hdi_2.5%\"]) & (hdi_multi[\"y_obs\"] <= hdi_multi[\"hdi_97.5%\"]), \n",
        "                    '#2F5597', '#C00000')\n",
        "    ax.scatter(hdi_multi[\"x_obs\"], hdi_multi[\"y_obs\"],\n",
        "            c = colors,\n",
        "            zorder = 2)\n",
        "\n",
        "    # 设置图例的颜色、形状、名称\n",
        "    legend_color = ['#2F5597', '#C00000']\n",
        "    handles = [plt.Line2D([0], [0], \n",
        "                        marker='o', \n",
        "                        color='w', \n",
        "                        markerfacecolor=color, markersize=10) for color in legend_color]\n",
        "    handles += [HDI]\n",
        "    labels = ['Within HDI', 'Outside HDI','95% HDI']\n",
        "\n",
        "    # 设置坐标轴名称、标题\n",
        "    ax.set_xlabel('Stress', fontsize=14) \n",
        "    ax.set_ylabel('Self control', fontsize=14)\n",
        "    fig.legend(handles=handles, labels=labels, loc='outside upper right')\n",
        "    fig.suptitle('Posterior Predictive Interval', fontsize=16)\n",
        "    sns.despine()\n",
        "\n",
        "# 绘制\n",
        "plot_ppi(model3_ppc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "9278EBAC74274EBBA3F603614FB7DF9B",
        "id": "31FC4F75AA364516B2D3355A27D39F8C",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 多元线性回归：增加交互项  \n",
        "\n",
        "* **交互作用：** 存在两个或两个以上自变量时，其中，自变量$X_1$对因变量$Y$的影响在自变量$X_2$的不同水平上有所差异  \n",
        "    \n",
        "    * 在本例中，我们可以假设，在不同的吸烟状况下，压力对自我控制的影响略有不同(体现在斜率上)  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s49oqb24we.png?imageView2/0/w/960/h/960)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "B435520B165044F8AE7D2E279B03AEAC",
        "id": "8B1CA2D03645490A9D4C9FB85727BCF3",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "那么我们可以加入交互项：$\\mu = \\beta_0 + \\beta_1 X_{1} + \\beta_2 X_{2} + \\beta_3 X_{1}X_{2}$  \n",
        "\n",
        "$$  \n",
        "\\begin{equation}  \n",
        "\\begin{array}{lcrl}  \n",
        "\\text{data:} & \\hspace{.05in} & Y_i|\\beta_0,\\beta_1,\\beta_2,\\beta_3,\\sigma & \\stackrel{ind}{\\sim} N(\\mu_i, \\; \\sigma^2) \\;\\; \\text{ with } \\; \\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2} + \\beta_3X_{i1}X_{i2} \\\\  \n",
        "\\text{priors:} & & \\beta_{0c}  &  \\sim N\\left(50, 10^2 \\right)  \\\\  \n",
        "                     & & \\beta_1  & \\sim N\\left(0, 10^2 \\right) \\\\  \n",
        "                     & & \\beta_2  & \\sim N\\left(0, 10^2 \\right) \\\\  \n",
        "                     & & \\beta_3  & \\sim N\\left(0, 10^2 \\right) \\\\  \n",
        "                     & & \\sigma & \\sim \\text{Exp}(0.6)  .\\\\  \n",
        "\\end{array}  \n",
        "\\end{equation}  \n",
        "$$  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "864CD8918B2F470EACAD5D93DBE803ED",
        "id": "C0F0EB0AF24F4447B3D2948088E79E67",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**在这个例子中，各参数的意义：**  \n",
        "\n",
        "* 当$X_{i2} = 0$时，$\\mu_i = \\beta_0 + \\beta_1X_{i1}$  \n",
        "\n",
        "    * 表示不吸烟情况下，自我控制分数随压力分数变化的情况，二者的关系可以被简化为一条直线。  \n",
        "\n",
        "* 当$X_{i2} = 1$时，$\\mu_i = \\beta_0 + \\beta_1X_{i1} + \\beta_2X_{i2} + \\beta_3X_{i1} = (\\beta_0 +\\beta_2) + (\\beta_1 + \\beta_3)X_{i1}$  \n",
        "\n",
        "    * 表示吸烟情况下，自我控制分数随压力分数变化的情况  \n",
        "    * 注意截距项和斜率项的变化  \n",
        "    * 此时压力对自我控制的影响为$(\\beta_1 + \\beta_3)$，体现了吸烟对这一关系的影响  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "7F23E4BBD0F64C0392856156F40D4A81",
        "id": "054D040027CC4E14B4BDB814152B4472",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**直接运行模型定义、MCMC采样、后验预测部分**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "_id": "28A556BF3D3749A6B601EAF15106DE08",
        "collapsed": false,
        "id": "D37BFCCCD1F940379F1E2DD92577E42F",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model() as model4:\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=50, sigma=10)        #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)         #定义beta_1\n",
        "    beta_2 = pm.Normal(\"beta_2\", mu=0, sigma=10)         #定义beta_2\n",
        "    beta_3 = pm.Normal(\"beta_3\", mu=0, sigma=5)          #定义beta_3\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                 #定义sigma\n",
        "\n",
        "    stress = pm.MutableData(\"stress\",df.stress, dims=\"obs_id\")      #stress是自变量压力水平\n",
        "    smoke = pm.MutableData(\"smoke\",df.smoke, dims=\"obs_id\")         #smoke是自变量吸烟水平\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + \n",
        "                                beta_1*stress + \n",
        "                                beta_2*smoke +\n",
        "                                beta_3*stress*smoke, dims=\"obs_id\")      #定义mu，将自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol, dims=\"obs_id\") "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "_id": "633BB97B162A4F4389DD1FBDC06B5ACD",
        "collapsed": false,
        "id": "412AB4627979404E81312C8663C059B3",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with model4:\n",
        "    model4_trace = pm.sample(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "                             tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                             chains=4,                     # 链数\n",
        "                             discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                             random_seed=84735)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "_id": "93F6850326BE4554B02837FDC9503963",
        "collapsed": false,
        "id": "CE18AABA06AC46B8B09A06FC3457DBEB",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with model4:\n",
        "    model4_ppc = pm.sample_posterior_predictive(model4_trace, random_seed=84735)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "6E96BD0CFE304B7BB529EC83E6824059",
        "id": "FB03BD47975A415C9982A9B78E180960",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 可视化后验预测结果  \n",
        "\n",
        "#### 后验预测的线性模型  \n",
        "\n",
        "我们可以分别查看不同吸烟情况下，压力对自我控制的预测情况  \n",
        "\n",
        "1. 当$X_{i2} = 0$时，$\\mu_i = \\beta_0 + \\beta_1X_{i1}$：不吸烟情况下，自我控制分数随压力分数变化的情况  \n",
        "\n",
        "2. 当$X_{i2} = 1$时，$\\mu_i = (\\beta_0 +  \\beta_2) + (\\beta_1 + \\beta_3)X_{i1}$：吸烟情况下，自我控制分数随压力分数变化的情况  \n",
        "\n",
        "* 根据吸烟情况($X_{i2} = 0/1$)，将后验参数的预测值分成两部分  \n",
        "\n",
        "* 代入公式，$\\mu = \\beta_0 + \\beta_1 X_{i1} + \\beta_2 X_{i2} + \\beta_3 X_{i1}X_{i2}$，获得对应的后验预测均值  \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "metadata": {
        "_id": "829864AAE2574A4A8AC6C8395610087D",
        "collapsed": false,
        "id": "C36BBAF702D74174ABADC127BD01DA5E",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 分别获取smoke=0及smoke=1时，压力变量对应的索引值\n",
        "no_smoke_index = xr.DataArray(df.query(\"smoke == 0\")[\"stress\"])[\"index\"].values\n",
        "smoke_index = xr.DataArray(df.query(\"smoke == 1\")[\"stress\"])[\"index\"].values"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "_id": "5EA2B4F03626478F9FDA544C190E594D",
        "collapsed": false,
        "id": "63B48C683B28405DBB304C2237ABBB2A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 根据上一步获得的索引值，将参数的后验预测分布分成两部分\n",
        "# 一部分为吸烟情况下，参数的后验预测分布、压力分数、自我控制分数\n",
        "# 另一部分为不吸烟情况下，参数的后验预测分布、压力分数、自我控制分数\n",
        "in_nosmoke_para = model4_trace.sel(obs_id=no_smoke_index)\n",
        "in_smoke_para = model4_trace.sel(obs_id=smoke_index)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "BD3E56DB4D274A5390C0CDAAB558F6B1",
        "id": "E93152923C364E9785475FA04BCC020D",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "根据吸烟情况将后验预测数据分成两份  \n",
        "\n",
        "* 查看 `in_nosmoke_para`:`constant_data`中的`smoke`值全为0，  \n",
        "* 相对的，`in_smoke_para`:`constant_data`中的`smoke`值全为1，"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "_id": "1D398840095C47909B8BB198E3F68E45",
        "collapsed": false,
        "id": "3AFF743A49904E04B9F558D9CE809B4A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "in_nosmoke_para"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 30,
      "metadata": {
        "_id": "D557C89308004325A9B86A847395D3B6",
        "collapsed": false,
        "id": "ED61A1D2373F4339B18CC1F2D8C2256E",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "in_smoke_para"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "E1BEFA9C33CB4E4499AC9197E0B3938A",
        "id": "9D279BF574984EBB8116D38B6684AC3E",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "根据回归公式计算不同吸烟条件下的回归预测模型  \n",
        "\n",
        "$$  \n",
        "\\mu = \\beta_0 + \\beta_1 X_{1} + \\beta_2 X_{2} + \\beta_3 X_{1}X_{2}  \n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "_id": "18289E7232854C968F4D8230FE250C9E",
        "collapsed": false,
        "id": "B68F66393B474D958B5CF02E101DD2D8",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 在不吸烟情况下，将参数与自变量进行组合，生成回归线\n",
        "in_nosmoke_para.posterior[\"y_model\"] = in_nosmoke_para.posterior[\"beta_0\"] +\\\n",
        "                                    in_nosmoke_para.posterior[\"beta_1\"] * in_nosmoke_para.constant_data[\"stress\"] +\\\n",
        "                                    in_nosmoke_para.posterior[\"beta_2\"] * in_nosmoke_para.constant_data[\"smoke\"] +\\\n",
        "                                    in_nosmoke_para.posterior[\"beta_3\"] * in_nosmoke_para.constant_data[\"smoke\"] * in_nosmoke_para.constant_data[\"stress\"]\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "_id": "7A92B82CBA7F4E38869722B5F9CA0077",
        "collapsed": false,
        "id": "139C057AEEB44D8995B3A2753E19F526",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 在吸烟情况下，将参数与自变量进行组合，生成回归线\n",
        "in_smoke_para.posterior[\"y_model\"] = in_smoke_para.posterior[\"beta_0\"] +\\\n",
        "                                    in_smoke_para.posterior[\"beta_1\"] * in_smoke_para.constant_data[\"stress\"] +\\\n",
        "                                    in_smoke_para.posterior[\"beta_2\"] * in_smoke_para.constant_data[\"smoke\"] +\\\n",
        "                                    in_smoke_para.posterior[\"beta_3\"] * in_smoke_para.constant_data[\"smoke\"] * in_smoke_para.constant_data[\"stress\"]\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {
        "_id": "872937B519894F67A20A78B9056CA805",
        "collapsed": false,
        "id": "0EE7F53F1ED84292A596DB1FC42F3836",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "plot_ppc_lm(in_nosmoke_para, in_smoke_para)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "0FD5780287B44F0999417ED27C83F890",
        "id": "B1B3FF681085498EA289F1FC6F69E2BA",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "我们可以将有交互项和无交互项两种模型生成的回归线进行对比：  \n",
        "\n",
        "- 我们可以发现，当不加入交互项时(model3)，在两种吸烟条件下，压力对于自我控制的预测是一致的，即两条深绿色线平行。  \n",
        "- 当加入交互项时(model4)，两种吸烟条件下，两条深绿色线不再平行。  \n",
        "- 然而，进一步的后验推断可以发现，交互项的效应并不可靠 (后验分布包含0)  \n",
        "\n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s4g6n11vv5.png?imageView2/0/w/640/h/640)  \n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "_id": "93E2FF7F03CF4770A9C318350BF6E03E",
        "collapsed": false,
        "id": "A9890D067D1542A39EABEFDA8860F313",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.plot_posterior(model4_trace, var_names = [\"beta_3\"])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {
        "_id": "D688CB09F9714ADDBA8122969554F0B7",
        "collapsed": false,
        "id": "8543EA53311242698D36978C3C49307A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.plot_forest(model4_trace, var_names = [\"beta_1\",\"beta_3\"], combined=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {
        "_id": "339C244135FD4BFA925903BE0A213927",
        "collapsed": false,
        "id": "B03D0533813B4072AB2D5C38860FD476",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.plot_forest(model4_trace, var_names = [\"beta_0\",\"beta_2\"], combined=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "1E40BDFCBC82446E86EB5612A4AD6FF0",
        "id": "D14590CC9AD04D6ABBD97454D8B0C4DF",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 后验预测分布  \n",
        "\n",
        "* 通过pymc，绘制当前模型对自我控制分数的后验预测分布"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {
        "_id": "8B01446558CE43E5827CBA9DFE2C2D94",
        "collapsed": false,
        "id": "1912C5CF3E3E43FEA57CCB7DA1759C01",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.plot_ppc(model4_ppc, num_pp_samples=500)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "21813C44E09D4F3089831C2A36CBD875",
        "id": "3DF91308AA8149CB93BBC0B5263773BE",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 评估后验预测  \n",
        "\n",
        "#### MAE  \n",
        "\n",
        "$$  \n",
        "\\text{MAE} = median\\left(|Y_i - \\hat{Y}_i|\\right)  \n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {
        "_id": "60DDF349608D4F608715D555B5CFEA4A",
        "collapsed": false,
        "id": "062EBB28BE214233B769F84B4573E248",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "model4_MAE = MAE(model4_ppc)\n",
        "print(f\"模型4(交互模型)MAE: {model4_MAE:.2f}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "E4DBE724AFF6425E82109862192D4F57",
        "id": "498F01E55CB14B4FB80E15A6C4A68D3D",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "#### 后验预测区间  \n",
        "\n",
        "* 使用`az.summary`，可以直接得到每个Y后验预测值95%的可信区间  \n",
        "* 将真实数据与其合并，判断其是否在对应的94%可信区间内，并进行绘图  \n",
        "\n",
        "> 注意：在真实数据中，一个X可能对应着多个Y，因此可信区间在绘制时会出现重叠"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "metadata": {
        "_id": "D1D6ACFA5E364D2EA4F5FEDDD73F486A",
        "collapsed": false,
        "id": "99C265EF45704B629025D1EAB3B08296",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "model4_outliers = counter_outlier(model4_ppc)\n",
        "print(f\"模型4(交互模型)超出95%hdi: {model4_outliers:.2f}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 40,
      "metadata": {
        "_id": "203B357341024A4E8FBD5E9F5FB760B0",
        "collapsed": false,
        "id": "30DC85604C1F412F8A96C475245CBD18",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# 绘制\n",
        "plot_ppi(model4_ppc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "308A177F40C4461289773F7188800F86",
        "id": "DC8B0A73AE3F41368594C60261CC41ED",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 总结  \n",
        "\n",
        "至此，我们已经讨论了四种模型，其汇总结果如下：  \n",
        "\n",
        "- 四个模型对于后验预测区间外观测值数量的预测得到了相同的结果，即均无法预测174个观测值中的9个。  \n",
        "- 而对于 MAE而言，以压力为预测变量的线性模型的MAE最低(3.41)，表明其预测性能最好。而以是否吸烟为预测变量的模型MAE最高，其预测性能最低(4.45)。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 41,
      "metadata": {
        "_id": "7AE9E42B16B24978B6C6ED2C8373C307",
        "collapsed": false,
        "id": "65ECCDB81D0748CEB6C8D0FC778C0313",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "summary_result = pd.DataFrame({\n",
        "    \"model_name\":[\"continuous model\",\"category model\",\"multi-linear model\",\"interaction model\"],\n",
        "     \"MAE\":[model1_MAE,model2_MAE,model3_MAE,model4_MAE],\n",
        "     \"HDI_num\":[model1_outliers,model2_outliers,model3_outliers,model4_outliers]\n",
        "     })\n",
        "\n",
        "summary_result "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "180597F207E8413BBCC07AD37ADD1CF6",
        "id": "442F1D34F8B44D6C9627C1E07E3319F8",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "🎯**开放性练习**  \n",
        "\n",
        "使用 `METU` 站点的数据进行线性模型分析。  \n",
        "\n",
        "1. 建立包含两个变量，以及交互项的回归模型  \n",
        "2. 使用 MAE 和 后验预测区间对模型进行评估。  \n",
        "4. 进行统计推断，解释模型中各回归系数代表的含义。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "76D99710B0A24F3EA9C8F86F215F3845",
        "id": "0DA141FDCB6847308C03F90CAADA22B1",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 模型评估与比较 (Model evaluation & comparison)  \n",
        "\n",
        "关于自我控制的主题，我们已经探索了多种研究假设和多个模型，你可能想知道哪个是最好的模型？  \n",
        "\n",
        "模型评估与比较(Model evaluation & comparison)的目的在于选择最好的模型。  \n",
        "\n",
        "什么需要模型比较？  \n",
        "1. 例如，比较 mode4 和 model3 可以帮助我们确定压力和吸烟之间的交互或调节作用。  \n",
        "2. 例如，比较 model3 和 model2(或model1) 可以衡量增加预测因子是否能提升模型的预测能力。  \n",
        "- 总之，模型比较的目的随着研究目的变化而变化。  \n",
        "\n",
        "|模型|参数|解释|  \n",
        "|-|-|-|  \n",
        "|model1|scontrol ~ stress|一般线性回归模型：自变量和因变量均为连续变量|  \n",
        "|model2|scontrol ~ smoke|一般线性回归模型：自变量为离散变量|  \n",
        "|model3|scontrol ~ stress + smoke|多元回归模型：自变量为连续变量和离散变量的组合|  \n",
        "|model4|scontrol ~ stress + smoke + stress:smoke|多元回归模型：自变量额外增加了两个自变量间的交互或调节作用|"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "39C8873E134A41D08E0B4D8DBA092FC2",
        "id": "D038B4D954BA45069E47501BE6C79D20",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 模型评估的绝对指标和相对指标  \n",
        "\n",
        "在上节课，我们讨论了三种评价模型的视角：  \n",
        "- 公平性(How fair)：模型在数据收集和分析的整个流程中的公正性。  \n",
        "- 错误程度(How wrong)：模型在**实践**中是否有效？即是否能够准确地预测**样本**数据。  \n",
        "- 准确性(How accurate)：模型是否反映**现实规律**？即是否能够准确地预测**样本外**数据。  \n",
        "\n",
        "而之前课程所讨论的评估标准(MAE和后验预测区间)都是**绝对指标**，即评估模型在**样本**数据上的预测能力。  \n",
        "\n",
        "而**相对指标**更多的是评估模型在**样本外**数据上的，这更有利于比较不同模型的预测能力。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "3CE9E67BF36D4776A2AA29FA949CE21B",
        "id": "3E70BEDF942D4C7C8568EF1C0F4D22B2",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**模型评估：对新数据的预测准确性**  \n",
        "\n",
        "为了评估模型对**样本外数据**的预测能力，我们有以下3种策略：  \n",
        "1. 通过新数据对模型进行评估。我们可以收集新的数据，并检验模型的预测能力。  \n",
        "2. 从已有样本中拿出一部分数据用来预测，这被称为**交叉验证(cross validation)** 方法。  \n",
        "3. 交叉验证的近似方法。例如，使用对交叉验证近似的**信息熵指标**来评估模型对样本外数据的预测能力。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "369E33BE8C5E442E92B86897732060E2",
        "id": "BD31649C0F364EEC8E8929C8A184EDA2",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "根据第一种策略收集新数据来检验模型的预测能力是一种理所当然的直觉。该策略包括两个主要步骤：  \n",
        "1. 训练模型：  \n",
        "   -  首先，我们根据现有的数据训练一个模型。如之前的例子，我们可以训练用压力和吸烟预测自我控制的线性模型。  \n",
        "2. 测试模型：  \n",
        "   -  在这阶段，我们需要收集新的数据，并检验模型的预测能力。例如，我们可以将新数据中测量的压力和吸烟数据带入模型，从而获得对于新数据的后验预测数据。最后，计算后验预测数据与新数据的MAE，从而评估模型对样本外数据的预测能力。  \n",
        "\n",
        "但心理学数据不同于其他学科的数据，它常受到**时间因素**和**抽样**的影响。  \n",
        "- 比如，个体的压力状态可能随着季节变化，因此在不同季节收集到的数据会受到时间的影响。  \n",
        "- 又例如，训练模型的数据来自理工科的学生，而测试模型的数据来自心理学的学生。这样抽样的差异也会影响模型对样本外数据的预测能力。  \n",
        "\n",
        "因此，一种更高效的方法是，一次性多收集一些数据，选择其中的一部分作为预测数据。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "D53E4FF71B9C400597A5A34B4752F6E0",
        "id": "B9AE664371394C46A0B14240D78C8B06",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 交叉验证(cross validation)  \n",
        "\n",
        "但问题在于，我们选择哪一部分数据作为预测数据呐？或者说，我们该如何有效的对数据进行抽取呐？  \n",
        "\n",
        "**交叉验证(cross validation)** 的目的就在于：提供不同的抽取预测数据的策略  \n",
        "- 其关键在于从已有样本中拿出一部分数据当作预测数据。  \n",
        "\n",
        "\n",
        "![](https://pic1.zhimg.com/80/v2-e9ad5ba61cda7ebd02848f336607eb70_1440w.webp)  \n",
        "\n",
        "\n",
        "> 资料来源：【绝对干货】机器学习模型训练全流程！- 知乎 https://zhuanlan.zhihu.com/p/184673895"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "A723C23C46674F1E8A9F6429A14EB668",
        "id": "FED7C4DD3C564577A71D99A6FB3F5271",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "常见的交叉验证策略：  \n",
        "1. 分半交叉验证 (Split-half cross-validation)  \n",
        "\t- 分半交叉验证将观测数据对半分成两部分，分别在不同的数据集上拟合模型，并在另外一半数据集上验证模型，最后再对比不同的模型在两份数据集作为验证集时的预测准确度。  \n",
        "2. K 折交叉验证 (K-fold cross-validation)  \n",
        "\t- K 折交叉验证把数据分成 K 分，其中一份作为训练集（拟合模型，对参数进行估计），其余的 K-1 分数据集作为验证集，总共重复这个流程 K 次。以 K 次验证结果的均值作为验证标准。  \n",
        "3. 留一法交叉验证 (Leave-one-out cross-validation)  \n",
        "\t- 留一法交叉验证是 K 折交叉验证的一个特例，当分折的数量等于数据的数量时，K 折留一法便成了留一法交叉验证。留一法交叉验证相较于普通的交叉验证方法，几乎使用了所有数据去训练模型，因此留一法交叉验证的训练模型时的**偏差 (bias) 更小、更鲁棒**，但是又因为验证集只有一个数据点，验证模型的时候**留一法交叉验证的方差 (Variance) 也会更大**。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "86A1642125D94E04BF66DE109850A976",
        "id": "F32A62A7BB3B43488EBFC38CCBC7F583",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**K 折交叉验证 (K-fold cross-validation)**  \n",
        "\n",
        "K 折交叉验证在分半交叉验证的基础上，将数据集分成 K 份(称为 CV-K)，其中一份作为测试集，其余 K-1 份作为训练集，重复这个流程 K 次。  \n",
        "\n",
        "K 折交叉验证，以 K 次测试结果的**均值**作为验证标准。例如，在压力-自我控制的例子中：  \n",
        "- 我们可以使用 K=5 折的交叉验证，将数据集分成 5 份，每次使用 5 份数据作为训练集，1份数据作为测试集。  \n",
        "- 对每一次迭代，我们使用 5 份数据训练模型，然后使用剩下的一份数据进行测试，并计算相应的MAE。  \n",
        "- 重复这个流程 5 次，然后取每次MAE测试结果的均值作为最终的测试结果。  \n",
        "\n",
        "![](https://pic3.zhimg.com/80/v2-ff846ee7eefdcd425e123d9d31b4d58a_1440w.webp)  \n",
        "\n",
        "> 资料来源：【绝对干货】机器学习模型训练全流程！- 知乎 https://zhuanlan.zhihu.com/p/184673895"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "44A080A89B104D498E12C9D22A831861",
        "id": "489B6E65F2864C3ABA5CF8C05A63DAF9",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**留一法交叉验证 (Leave-one-out cross-validation)**  \n",
        "\n",
        "留一法交叉验证是 K 折交叉验证的一个特例，当分折的数量K等于数据的数量n时，K 折留一法便成了留一法交叉验证。  \n",
        "- 留一法交叉验证相较于普通的交叉验证方法，几乎使用了所有数据去训练模型。  \n",
        "- 留一法交叉验证 (Leave-one-out cross-validation)的缩写为 loo-cv，或者 loo。  \n",
        "\n",
        "![](https://www.baeldung.com/wp-content/uploads/sites/4/2022/05/loso.png)  \n",
        "\n",
        "> 资料来源：https://www.baeldung.com/cs/cross-validation-k-fold-loo"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "661077D68096433EAEAF059085D90742",
        "id": "659D185F93DA49BD9122F43E5B014F3A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### ELPD (Expected log-predictive density)  \n",
        "\n",
        "留一法交叉验证 LOO (包括之前的交叉验证方法)是用于评估模型在**未知数据**上预测能力的思想框架，其本身并不提供具体的统计指标。  \n",
        "\n",
        "**ELPD** (Expected log-predictive density) 是 LOO 方法的具体实现，以对数似然函数作为统计指标。  \n",
        "\n",
        "其计算步骤：  \n",
        "- 同 K 折交叉验证一样，首先将数据集分成 n 份，n为数据总的数量。  \n",
        "- 利用 n-1 份数据去训练模型，得到后验模型 $p(\\theta_{-i}|y_{-i})$。  \n",
        "- 使用剩下的一份数据作为测试数据 $y_{i}$，计算后验预测模型 $p(y_{i}|y_{-i})$。  \n",
        "- 重以上过程，重复 n 次，得到 n 个后验预测模型,并计算其对数化后的期望值 $E(log(p(y_{i}|y_{-i})))$。  \n",
        "\n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s4bmgg7han.png?imageView2/0/w/640/h/640)  \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "276B0263F1FC4E319259A46AE6BA1141",
        "id": "D320603E9DE84AFB99B8E132825855CB",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "在实际操作中，我们通过 `ArViz` 的函数`az.loo`计算 $ELPD_{LOO-CV}$。  \n",
        "\n",
        "* 在`az.loo`返回的值中，`elpd_loo` 为$E(log(p(y_{i}|y_{-i})))$，  \n",
        "* `elpd_loo`越高表示模型的预测值越精确  \n",
        "\n",
        "注意：由于 $ELPD_{LOO-CV}$ 的计算量也比较大，ArViz 会使用 Pareto Smooth Importance Sampling Leave Once Out Cross Validation (PSIS-LOO-CV) 来近似 $ELPD_{LOO-CV}$。  \n",
        "\n",
        "PSIS-LOO-CV 有两大优势：  \n",
        "1. 计算速度快，且结果稳健  \n",
        "2. 提供了丰富的模型诊断指标  \n",
        "\n",
        "注意：  \n",
        "- 要计算 elpd_loo 需要在采样 `pm.sample` 中加入 `idata_kwargs={\"log_likelihood\": True}`  \n",
        "- 或者，在模型采样完成后计算对数似然，即 `with model:  pm.compute_log_likelihood(model_trace)`。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 42,
      "metadata": {
        "_id": "9A0AD2AD255246D3A25BE05C8FA49F2B",
        "collapsed": false,
        "id": "99E272BB5FA54A1CA5A6071025AA09FD",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "with pm.Model() as model1:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=45, sigma=10)        #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)         #定义beta_1\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                 #定义sigma\n",
        "\n",
        "    x = pm.MutableData(\"smoke\",df.stress)                     #x是自变量压力水平\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + beta_1*x)            #定义mu，讲自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol)\n",
        "\n",
        "\n",
        "with pm.Model() as model2:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=45, sigma=10)        #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)         #定义beta_1\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                 #定义sigma\n",
        "\n",
        "    smoke = pm.MutableData(\"smoke\",df.smoke)            #x为是否吸烟，0 1编码。\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + beta_1*smoke)  #定义mu，讲自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol)\n",
        "with pm.Model() as model3:\n",
        "\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=50, sigma=10)          #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)           #定义beta_1\n",
        "    beta_2 = pm.Normal(\"beta_2\", mu=0, sigma=10)           #定义beta_2\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                   #定义sigma\n",
        "\n",
        "    stress = pm.MutableData(\"stress\",df.stress)                           #自变量1:压力水平\n",
        "    smoke = pm.MutableData(\"smoke\",df.smoke)                              #自变量2:是否吸烟\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + beta_1*stress + beta_2*smoke)    #定义mu，自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol) \n",
        "with pm.Model() as model4:\n",
        "    beta_0 = pm.Normal(\"beta_0\", mu=50, sigma=10)        #定义beta_0          \n",
        "    beta_1 = pm.Normal(\"beta_1\", mu=0, sigma=10)         #定义beta_1\n",
        "    beta_2 = pm.Normal(\"beta_2\", mu=0, sigma=10)         #定义beta_2\n",
        "    beta_3 = pm.Normal(\"beta_3\", mu=0, sigma=5)          #定义beta_3\n",
        "    sigma = pm.Exponential(\"sigma\", 0.6)                 #定义sigma\n",
        "\n",
        "    stress = pm.MutableData(\"stress\",df.stress)      #stress是自变量压力水平\n",
        "    smoke = pm.MutableData(\"smoke\",df.smoke)         #smoke是自变量吸烟水平\n",
        "    mu = pm.Deterministic(\"mu\", beta_0 + \n",
        "                                beta_1*stress + \n",
        "                                beta_2*smoke +\n",
        "                                beta_3*stress*smoke)      #定义mu，将自变量与先验结合\n",
        "\n",
        "    likelihood = pm.Normal(\"y_est\", mu=mu, sigma=sigma, observed=df.scontrol) \n",
        "\n",
        "#===========================\n",
        "#     注意！！！以下代码可能需要运行15分钟左右\n",
        "#===========================\n",
        "with model1:\n",
        "    model1_trace = pm.sample(draws=5000,                  # 使用mcmc方法进行采样，draws为采样次数\n",
        "                      tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                      chains=4,                     # 链数\n",
        "                      discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                      idata_kwargs={\"log_likelihood\": True},\n",
        "                      random_seed=84735)\n",
        "with model2:\n",
        "    model2_trace = pm.sample(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "                        tune=1000,                   # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                        chains=4,                    # 链数\n",
        "                        discard_tuned_samples= True, # tune的结果将在采样结束后被丢弃\n",
        "                        idata_kwargs={\"log_likelihood\": True},\n",
        "                        random_seed=84735)  \n",
        "with model3:\n",
        "    model3_trace = pm.sample(draws=5000,            # 使用mcmc方法进行采样，draws为采样次数\n",
        "                      tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                      chains=4,                     # 链数\n",
        "                      discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                      idata_kwargs={\"log_likelihood\": True},\n",
        "                      random_seed=84735)\n",
        "with model4:\n",
        "    model4_trace = pm.sample(draws=5000,                   # 使用mcmc方法进行采样，draws为采样次数\n",
        "                                tune=1000,                    # tune为调整采样策略的次数，可以决定这些结果是否要被保留\n",
        "                                chains=4,                     # 链数\n",
        "                                discard_tuned_samples= True,  # tune的结果将在采样结束后被丢弃\n",
        "                                idata_kwargs={\"log_likelihood\": True},\n",
        "                                random_seed=84735)    "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 43,
      "metadata": {
        "_id": "3EB067CA651442A78F7346CACAC53310",
        "collapsed": false,
        "id": "8E3359468F2B438B941D1717E3016E0A",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "az.loo(model4_trace)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "1E007CED1FE14C5890758363889A6FC3",
        "id": "BF5AF5A67F1D49FD9D02CE233BE819DE",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "模型`elpd_loo`的结果为-545.62  \n",
        "\n",
        "- 然而，仅凭单个值，并不能反映模型的预测精确程度。  \n",
        "- 虽然 **ELPDs 无法为任何单一模型的后验预测准确性提供可解释的度量，但它在比较多个模型的后验预测准确性时非常有用**。  \n",
        "\n",
        "我们可以通过 `arviz.compare` 方法来对比多个模型的 elpd。从下面结果可见：  \n",
        "- 模型1的 elpd_loo 最大，表明它对**样本外数据**的预测性能最好。  \n",
        "- 而模型2的 elpd_loo 最小，表明它的预测性能最差。  \n",
        "- 并且这些结果与我们通过 MAE 和 后验预测区间得到的判断一致。  \n",
        "\n",
        "需要注意的是：  \n",
        "- arviz 提供的结果包括了 elpd se，这使得我们可以判断两个模型的预测差异 elpd_diff 是否超过两至三个标准误se。  \n",
        "- 从现在的结果看，模型1，3，4的差异很小；而模型2的和其他模型的差异 (19.15) 已经超过了两个标准误se。表明该模型的预测性能显著的低于其他三个模型。"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 44,
      "metadata": {
        "_id": "1C9EBCCA426A4C8294CDD7ED0D40E11F",
        "collapsed": false,
        "id": "6BCD16F940454059BDB0A8CB98ED9E0C",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "comparison_list = {\n",
        "    \"model1(contiunous)\":model1_trace,\n",
        "    \"model2(category)\":model2_trace,\n",
        "    \"model3(multivariate )\":model3_trace,\n",
        "    \"model4(interaction)\":model4_trace,\n",
        "}\n",
        "az.compare(comparison_list)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 45,
      "metadata": {
        "_id": "CD740120DBFA47AC9548AED747B8369F",
        "collapsed": false,
        "id": "D1792B9CACB744E2A63C01B9D1038420",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": [],
        "trusted": true
      },
      "outputs": [],
      "source": [
        "summary_result"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "17BA4AFC1FCF4C1495BA651B224862F1",
        "id": "34954DC93C354D37BB773DE2683988FD",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "🤔思考：  \n",
        "\n",
        "1. 如果你的目标是在不控制任何其他因素的情况下探索吸烟与自我控制之间的关系，你会使用哪种模型？  \n",
        "    \n",
        "    - 考虑到模型3(多元回归模型)优于模型2(仅包含吸烟)--因此，选择模型3可能更好地反映吸烟的真实效应，而不是忽略压力的主要作用。  \n",
        "        \n",
        "2. 如果你的目标是最大限度地提高模型的预测能力，而在模型中只能选择一个预测因子，您会选择压力还是是否吸烟？  \n",
        "    \n",
        "    - 由于模型1优于模型2，如果仅选择一个预测变量的话，选择压力能获得对于自我控制更好的预测。  \n",
        "        \n",
        "3. 这四个模型中，哪个模型的**总体预测结果**最好？  \n",
        "    - 模型1以微弱优势超过了使用所有预测因子的模型4。这表明，在建立模型的过程中，预测因子并不是越多越好。  \n",
        "    - 事实上，模型4比模型3还差一点，这表明压力和吸烟之间的交互效应很弱，加入两者的交互项，会减弱模型的预测能力。  \n",
        " \n",
        " 因此，为了简单高效，我们更有理由选择模型1或者模型3。  \n",
        " \n",
        " ![](https://th.bing.com/th/id/R.31c49d0bc73e477eda2cf52fef1b859f?rik=fdM9GvPWnNMg4Q&riu=http%3a%2f%2fwww.esafety.cn%2fblog%2fUploadFiles%2f2018-7%2f51155444445.jpg&ehk=O2MXHgnsmdvS68QMhw6CaIw82aHg2E0q%2fSKKtutAZjk%3d&risl=&pid=ImgRaw&r=0)  \n",
        " \n",
        " > 资料来源: http://www.esafety.cn/Blog/u/9490/archives/2018/154367.html"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "17368D2259514152A39C7E07EB11D00B",
        "id": "940DA9B666BD464EB7B98CD305D04FA9",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**如何增加后验预测的能力？**  \n",
        "\n",
        "往往有两种通用的策略：  \n",
        "- 收集更多的数据：  \n",
        "    - 虽然收集更多的数据并不能改变 Y 和 X 之间潜在的微弱关系，但它可以改善模型对参数后验的估计，从而改进模型的后验预测。  \n",
        "    - 一些可以考虑的手段，增加数据量(n)，扩大x的范围，扩大后验预测的范围。  \n",
        "- 加入更多的预测变量：  \n",
        "    - 例如，在自我控制能力方面，压力比吸烟有更好的预测。  \n",
        "    - 此外，我们可以通过加入 2 个或更多的预测因子来改进我们的模型及其对 Y 的预测（即同时将压力和吸烟作为预测因子加入模型）。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "D80EE8DC5B9E4DBB8B272571EBB29DA4",
        "id": "D68B37FE67FD47D68221A3C9AF3C72A1",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 偏差方差权衡(bias-variance trade-off)  \n",
        "\n",
        "模型比较的例子和结果说明：  \n",
        "- 如果我们的目标是建立一个能准确预测响应变量 Y 的模型，那么我们就需要包含足够多的预测因子，这样我们就能获得有关 Y 的充分信息。  \n",
        "- 然而，加入太多的预测因子，有可能适得其反。  \n",
        "\n",
        "偏差方差权衡就反映了这样一个事实：  \n",
        "- 当模型越复杂，可能对于样本数据有更好的解释，然而无法解释样本外的数据。  \n",
        "- 当模型太简单时，模型可能照顾到不同的样本，但对于所有样本的解释力都一般。  \n",
        "\n",
        "![Image Name](https://vitalflux.com/wp-content/uploads/2020/12/overfitting-and-underfitting-wrt-model-error-vs-complexity-768x443.png)  \n",
        "\n",
        "资料来源：https://vitalflux.com/overfitting-underfitting-concepts-interview-questions/"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "1D4283AD689843FAAC80213F9A26778C",
        "id": "2E42D42214AD49C8A60DAE1582F0C996",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "模型评估的核心在于模型捕捉到了数据中的关键模式，既非太简单而错过数据中有价值的信息(**欠拟合, underfitting**)，也不会太复杂从而将数据中的噪音加入到模型中(**过拟合, overfitting**)。  \n",
        "\n",
        "**欠拟合(underfitting)**  \n",
        "\n",
        "* 欠拟合的模型在当前样本的数据拟合效果不好，且其泛化能力(模型在当前样本外新的数据上的预测的准确度)也同样不佳。  \n",
        "* 导致欠拟合的原因  \n",
        "  * 数据特征较少  \n",
        "    * 数据特征指的是数据的属性，比如第一部分中展示的数据的各个变量就是数据的特征。在所有变量都能独立地对目标变量做出解释的前提下，数据特征越多，数据拟合程度越好。  \n",
        "  * 模型复杂度过低  \n",
        "    * 模型的复杂度代表模型能够描述的所有函数，比如线性回归最多能表示所有的线性函数。  \n",
        "    * 模型的复杂度和模型的参数数量有关，一般来说，模型参数越多，复杂度越高，模型参数越少，复杂度越低。  \n",
        "\n",
        "**过拟合(overfitting)**  \n",
        "\n",
        "* 模型在当前样本的数据上的拟合程度极好，但是泛化能力也较差。  \n",
        "* 模型把训练样本学习地“太好了”，把样本自身地一些噪音也当作了所有潜在样本都会具有的一些性质，这样就会导致其泛化性能下降。  \n",
        "* 导致过拟合的原因  \n",
        "  * 当前样本的噪音过大，模型将噪音当作数据本身的特征  \n",
        "  * 当数据的有些特征与目标变量无关，这些特征就是噪音，但它也可能被误当作数据特征，这就会造成模型过拟合  \n",
        "  - 样本选取有误，样本不能代表整体  \n",
        "  - 模型参数太多，模型复杂度太高  \n",
        "\n",
        "\n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/s4bp7piw04.png?imageView2/0/w/640/h/640)  \n",
        "\n",
        "\n",
        "资料来源：https://blog.csdn.net/weixin_43378396/article/details/90707493"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "C92F7C50883E465F9FE447381B96F0D0",
        "id": "5511BB7FABA84B7F8DFA1927585B238B",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "### 如何避免欠拟合  \n",
        "\n",
        "- 增加数据的特征  \n",
        "\n",
        "- 增加模型复杂度  \n",
        "\n",
        "### 如何避免过拟合  \n",
        "\n",
        "- 选择更具代表性的数据  \n",
        "\n",
        "- 降低模型复杂度  \n",
        "\n",
        "\n",
        "**问题的本质在于：模型与数据真实的生成模型匹配**  \n",
        "\n",
        "为了选择一个能够在过拟合和欠拟合之间的达到平衡的最佳模型，就需要进行模型评估、比较和选择。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "11121B94841F40B4AD09F0DE3570342F",
        "id": "2EA18F9E553D406B94265819EE6F2A36",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## workflow  \n",
        "\n",
        "**概念区别：模型诊断(model diagnostic) vs. 模型验证(model validation)**  \n",
        "\n",
        "在贝叶斯的workflow中，  \n",
        "\n",
        "- MCMC评估是对MCMC进行检验，在Gelman et al (2020)中也被称为计算过程的验证(validate computation)，其核心在于确定MCMC算法采样得到的样本是否足以提供目标分布的精确近似。  \n",
        "\n",
        "- 模型评估则是指对模型是否公平性、有效性、可信性进行评估，既可以是对单个模型，也可以是对多个模型进行。  \n",
        "\n",
        "理想情况下，两种评估结果都是良好的。然而，情况并非总是如此。  \n",
        "  - 我们可能有一个很好的模型框架，但是MCMC链的结果却不稳定，进而导致不准确的后验近似。  \n",
        "  - 或者刚好相反，我们可以获得很好地后验近似，但模型本身却并不反映任何事实。  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/image/rkz1dqyo2e.png?imageView2/0/w/960/h/960)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "6BB7D09582954DA4B4C70FD2D705BD40",
        "id": "A7F3DAC1EB934D2D9DE27D2717F148B3",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "🎯**开放性练习**  \n",
        "\n",
        "请使用 `METU` 站点的数据完成整个贝叶斯 workflow 的分析。  \n",
        "1. 根据研究问题或假设建立多个回归模型。需要自行决定每个回归模型的预测变量和先验。  \n",
        "2. 使用 pymc 或者 bambi 拟合模型。  \n",
        "3. 进行模型评估，包括通过 LOO 方法进行模型比较，以及对最优模型进行后验预测检验和后验预测数据的分析。  \n",
        "4. 进行统计推断，报告参数的后验估计结果(参数后验的均值以及95HDI)。"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "3CF7E4FB7F0B48F79EB864DAC5DA5314",
        "id": "48B2D6B27347407F911DE5118A3B6BD0",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "## 总结  \n",
        "\n",
        "本节课以不同模型的评估为视角，介绍了模型评估和比较的思想  \n",
        "\n",
        "至此我们也对贝叶斯分析的全流程(Bayesian workflow)有了一个初步的认识  \n",
        "\n",
        "之后的课程中，我们会不断实践这个流程，让大家感受贝叶斯分析的魅力。  \n",
        "\n",
        "![Image Name](https://cdn.kesci.com/upload/image/rkz1ehen1l.png?imageView2/0/w/960/h/960)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_id": "8E9B15CA23FB4440B064FC86D9D553C0",
        "id": "322B8D87F92545048F0A218905CBE5E1",
        "jupyter": {},
        "notebookId": "655b69f19f41eb29b96aa415",
        "runtime": {
          "execution_status": null,
          "is_visible": false,
          "status": "default"
        },
        "scrolled": false,
        "slideshow": {
          "slide_type": "slide"
        },
        "tags": []
      },
      "source": [
        "**补充：**  \n",
        "\n",
        "之前讨论过模型评估中两种类型的指标：  \n",
        "- 绝对指标，衡量模型对于样本的预测能力。  \n",
        "- 相对指标，衡量模型对于样本外数据的预测能力，也考虑了模型的复杂度。  \n",
        "\n",
        "这两种指标包含了多种具体的统计值：  \n",
        "\n",
        "模型拟合优度的方法包括：  \n",
        "- MAE or MSE(mean square error)  \n",
        "- 对数似然 (log likelihood)  \n",
        "- $R^2$  \n",
        "\n",
        "模型预测进度的方法包括：  \n",
        "- AIC  \n",
        "- DIC  \n",
        "- WAIC  \n",
        "- LOO-CV  \n",
        "\n",
        "\n",
        "|                    | AIC                                  | DIC                                      | WAIC       | LOOCV           | BIC                                  |  \n",
        "| ------------------ | ------------------------------------ | ---------------------------------------- | ---------- | --------------- | ------------------------------------ |  \n",
        "| 适用框架           | 频率论                               | 贝叶斯                                   | 贝叶斯     | 贝叶斯          | 贝叶斯/频率论                        |  \n",
        "| 偏差（deviance）   | 最大似然参数 $\\theta_mle$ 的对数似然 | 贝叶斯参数均值 $\\bar{\\theta}$ 的对数似然 | LPPD       | $ELPD_{LOO-CV}$ | 最大似然参数 $\\theta_mle$ 的对数似然 |  \n",
        "| 矫正（correction） | 参数数量                             | 似然的变异                               | 似然的变异 |     由于采用 LOO-CV 思想，因此不需要矫正            | 参数数量+数据数量                    |"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.5.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}
